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Abstract 

We formulate dynamical rate equations for physical processes driven by a combination of diffu- 
sive growth, size fragmentation and fragment coagulation. Initially, we consider processes where 
coagulation is absent. In this case we solve the rate equation exactly leading to size distributions of 
Bessel type which fall off as exp(— x'^/^) for large x-values. Moreover, we provide explicit formulas 
for the expansion coefficients in terms of Airy functions. Introducing the coagulation term, the 
full non-linear model is mapped exactly onto a Riccati equation that enables us to derive various 
asymptotic solutions for the distribution function. In particular, we find a standard exponential 
decay, exp(— x), for large x, and observe a crossover from the Bessel function for intermediate 
values of x. These findings are checked by numerical simulations and we find perfect agreement 
between the theoretical predictions and numerical results. 
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I. INTRODUCTION 

Since the pioneering work of Smoluchowski dating back to the beginning of the last 
century, the hterature on coagulation and fragmentation processes has grown considerably. 
Smoluchowski original coagulation equation [Jl^ provides a mean field description of clusters 



that coalesce by binary collisions with a constant rate. Scaling theory and exactly solvable 
models in the kinetics of irreversible aggregation has recently been reviewed in 3[. 

Fragmentation and coagulation was first considered as combined processes in P], and 



mean field type of coagulation-fragmentation models have subsequently been used in a di- 

cluster formation in 



verse range of applications, including polymer kinetics 0, aerosols 

astrophysics 3] and animal grouping in biology 0, 0] . We refer to |lo|, (and references 
therein) for a survey of the progress in the study of coagulation-fragmentation process. 

Recently, we have suggested a mean field model describing the dynamics of coherent 
structures, like ice crystals and structural elements of biomolecules, which grow or shrink 
randomly due to the diffusive motion of their boundaries and are subject to occasional 
fragmentation jl^ . llS^ . In this paper, we present an extension of the model to account for 
coagulation processes as well. The extension can be considered as a generalization of the 
Smoluchowski coagulation-fragmentation equation to include processes where size diffusion 
is important. We emphasize that our approach differs from the approach in e.g. where 
the clusters represent particles immersed in a gas or liquid type of medium and the diffusive 
term added to the coagulation-fragmentation equation represents the random movement of 
the center of mass of each cluster. In contradistinction, we are focusing on systems where 
diffusion operates in the size space rather than in real space. 

Diffusive growth processes are encountered in numerous physical systems. For instance, 
the process behind grain growth in ice or metallurgical systems can effectively be described 
as a size diffusive process, where the diffusion constant depends on the surface tension and 
the mobility of the grain boundaries (in the absence of extrinsic drag forces resulting from 
i. e. impurities in the material), see i. e. 12, 13, ll^ ll^ ll^- Size diffusion also appears 



in natural conjunction with convective growth in systems where the clusters in concern are 
coupled to a medium mediating the addition or subtraction of monomeric units. Crystal 
growth dynamics by reversible solute addition has for instance been considered in [19], in 
the limit of no coagulation and fragmentation. The combined process of fragmentation. 
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2l| to describe the abrupt transi- 



size diffusion and convection has been studied in 2^ 
tions from the growing to the shrinking state of microtubules. In both cases, the effective 
convective and diffusive growth can be directly related to the rate equations for monomer 



tions in polypeptide systems 



attachment /detachment from and to the medium (see i.e. ^^). Finally, structural transi 



22j | can be modeled by a continuous master equation with a 



size diffusive term representing the random fluctuations of the ordered-disordered interface 



121. 



Although we are not aware of previous work where diffusive growth appears in conjunc- 
tion with cluster coagulation, we believe -in light of the examples given above- that such 
extension to Smoluchowskis original fragmentation-coagulation model is natural and should 
find application in several physical systems. Our aim in this paper is therefore to elucidate 
the effect of size diffusion in simple models of coagulation-fragmentation processes. 

To proceed, let N{x,t) denote the number of clusters of size x at time t. The general 
form of the coagulation-fragmentation model then reads Q| 

mix, t) = [m]co^^ + [dtN]^^^^. (1) 

where 

1 r 



\dtN]cos.R = - Kix- x\ x')N{t, x - x')N{t, x')dx' - N{t, x) / K{x, x')N(t, x')dx' 
2 JO JO 



and 



[^<^]fraff = -N{x, t) J F{x - x', x)dx + 2 J F{x, x')N{x + x' , t)dx' . 



rrag 

The microscopic details of a given physical model is embedded in the kernels K{x, x') > and 
F{x, x') > (symmetric in x and x') which represent, respectively, the rate of coagulating 
two clusters of size x and x' into one cluster of size x + x' and the rate of fragmentation of a 
cluster of size x + x' into two clusters of size x and x'. The possible exchange of monomeric 
units with an external medium can be modeled by adding a convective and diffusive term 
to the equation which then reads, 

dtN{x, t) = -v{x)d^N{x, t) + D{x)dlN{x, t) + [9iiV]coag + [dtN]^.,^^. (2) 

Here, the first term describes the average drift of the cluster sizes due to monomeric exchange 
with the medium and the second term represents the random fluctuations superimposed on 
this drift (size diffusion) . In the following we will examine analytical and numerical aspects of 



this equation with the simple choice of size-independent parameters: p (x) = D, F{x, x') = f 
and K{x,x') = [3 and setting f = 0. The model presented in [ij, |l3| for the dynamics 
of ice crystals and structural building blocks in biomolecules, corresponds to case of no 
coagulation (3 = 0. In these cases the systems are considered closed, so the f = condition 
follows from excluding any external driving. When the system is in a solution allowing for 
the solute-solvent exchange of monomeric units (as in i.e. jiol. l2ol. 21 1) the f = condition 



corresponds to a situation where the solute (clusters) are in thermodynamical equilibrium 
with the solvent. The existence of steady states in the model with no coagulation has 
recently been identified for a wider class of fragmentation kernels j^], where connections to 



the pure fragmentation equation 



24| and the so-called "shattering" transition j25| related to 



the formation of dust particles, is further discussed. In the absence of diffusion, the model 
has been applied to the kinetics of reacting polymers in 
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27l |. where the uniqueness of 



solutions and convergence to equilibrium in the limit t ^ oo have been proved. As will 
be demonstrated, the inclusion of a diffusion term in the equation may lead to solutions 
radically different from the ones found in [2^ ■ 

Our paper is organized as follows. In section 2, we derive the exact solution to the model 
in the case of no coagulation (D > 0, / > 0, /3 = 0). Section 3 and 4 is devoted to 
the discussion of the expansion coefficients entering this solution. In particular, we give 
an explicit formula for the coefficients in terms of a orthogonal set of Airy functions and 
discuss their asymptotic behavior. A moment analysis of the solution in the large time limit is 
carried out in section 5. In section 6, we focus on the full diffusion-coagulation-fragmentation 
model. It is shown that the equation can be mapped to a Riccati equation, and we discuss 
its solutions in the case of no coagulation (/5 = 0, section 7) and no diffusion {D = 0, section 
8). Here, the solution to the diffusion- fragmentation model found in section 2-5 becomes 
useful to elucidate the structure of the solutions to the Riccati equation. In the section 
9 we return to the full diffusion-fragmentation-coagulation model. We demonstrate that 
the equation possesses a transition point for the coagulation term, where the distribution 
crosses over from the Bessel behavior found in the pure fragmentation-diffusion case to an 
exponential. In the final section, we discuss two numerical implementations of the model. 
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II. SOLUTION TO THE FRAGMENTATION-DIFFUSION EQUATION 

In the following three sections, we expand on the results of a recent letter on the 
diffusion-fragmentation equation. In particular, we shall carefully study the structure of 
the analytical solution of the diffusion-fragmentation equation. We shall derive orthogo- 
nality relations of the functions entering the solution and find the corresponding expansion 
coefficients. The general equation we wish to solve has the form 

dtN{x, t) = DdlN{x, t) - fxN{x, t) + 2f N{x', t)dx', 



which is a particular case of Eq. Q with size independent diffusion and fragmentation 
terms, and with no coagulation {(3 = 0). Below we shall consider the boundary conditions 
A^(0, t) = and N{x, t) — > for a; — >^ oo for all times t, including the initial condition at 
t = 0. 

We rescale time and the cluster size 

x^x/xo and t^t/(fxoy^, 

where xq = {D/ fY^^. A subsequent differentiation with respect to x leads to the equation: 

dtd^N{x, t) = dlN{x,t)- 3N{x, t) - xd^N{x, t) . (3) 

It is fairly easy to solve Eq. Q by separation of variables. Taking N{x,t) = T{t)X{x) we 
have 

f(t) = XT{t), X"'{x) - 3X{x) - xX\x) = XX'{x), (4) 

where A is a separation constant. We shall simplify the notation by shifting the variable x 
such that the separation constant gets included, x ^ x — X. The equation for X can be 
solved (with the boundary condition X{x) — > for x — > oo) by means of a Fourier transform, 
with the result 

X{x) = - dke e^^^^^ 1^ = B{x). (5) 

2 J— oo 

Although X is positive we need to consider X{x) for negative values of x in equation (jH)). 
By analytical continuation the function X{x) is well defined for x < 0. The function B is 
related to the Airy integral 

A[x) = - dk e*fc"+*'= /3 = / dk cos {kx + — ] , (6) 
2 J-oo Jo V 3 / 



which has been discussed eg. in reference j28|, see in particular pp. 188-190. In the panels 
of Fig. [T]we show B{x) for positive and negative values of x, respectively. 
By differentiating A{x) twice we get 



B(x) 



-dl A{x) 



(7) 



Collecting the results obtained above, and inserting once again the separation constant, the 
solution of Eq. ^ can now be written 



N{x,t) = Y. Cne^-'B{x + Xn] 



(8) 



In case the eigenvalues A„ are in a continuous range the sum should be replaced by an 
integral. We shall here impose an absorbing boundary conditions for small clusters, which 
implies that the probability of clusters with zero size vanishes at all times, i. e. 



Ar(0,t) = 0. 



(9) 



We implement this condition by requiring -B(A„) = 0. Thus the eigenvalues must be zeros 
of the function B{x). 

To find possible zeros of B{x) we need some properties of the functions A{x) and B{x). It 
turns out that A{x) can be expressed in terms of Bessel functions by use of a deformation 
of the contour in Eq. (jHl), with the result 



A{x) 



J- Ki for x > 

V 3 3 I 3 / 



3^/1x1 



Ji 

3 



J 1 



2\x 



3/2' 



for X < 0. 



(10) 



The absolute signs in the last line should be noticed. The function A{x) is positive for x > 
and oscillates for x < 0, where it has an infinite number of zeros. 

rn 

We can now differentiate A(x), and using well known functional relations i32l for the 

n n 

Bessel functions 1281 , one obtains the result 



A'{x) 



X / 2x3/2 \ 



vr 



-X 



J_2 

3 



'2|x|3/2^ 



J2 

3 



'2|xP/2^ 



for X < 0. 



(11) 
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FIG. 1: The function B(x) (x is dimensionless) and its intersections with y = 0. Note that all 
intersections appears for x < only. Inset: A zoom of the function for positive x. 

Performing another differentiation gives 

A"{x) = (^^^ for X > 0, 

'2|xp/2\ ^ /2|x|3/2^ 



TT 



h 13 j + '-i jj - < 0- (12) 

Comparison with fIlOp shows that A{x) satisfies the simple second order differential equation 
found by Stokes ^ 

dlA{x) = xA{x), (13) 

which will play a crucial role in the following. Comparing this equation to (|7j) we see the 
remarkably simple relation between B{x) and A{x), 

B{x) = -xA{x). (14) 

From well known properties of the Bessel functions K and J we see that the function 
A{x) has zeros for x = A„ < with 



0. 



(15) 



These zeros can be computed easily numerically. The first few are given approximately by 
Ai = -2.338, A2 = -4.088, A3 = -5.521, A4 = -6.787, A5 = -7.945, Ag = -9.023. (16) 
It should be mentioned that from the asymptotic form of the Bessel functions the zeros are 
approximately given by 



A. 



37r/ 3M2/3 



, n + 
2 V 4 



(17) 



This expression is valid for large n. However, the expression provides a surprisingly accurate 
estimate of A„+i even for small n's. For example, for = one finds that Eq. ()17|) gives 
—2.32, as compared to the numerically obtained more accurate value Ai = —2.338. For 
n = 3 one obtains ~ —6.785 to be compared to the more accurate A4 = —6.787. 

We now return to the boundary condition Q which is to be implemented on the solution 
(jH)). This requires 

00 

iV(0,t) = ^ C„5(A„) = 0. (18) 

n=0 

Taking into account Eq. (fT^ giving the relation between A{x) and B{x) as well as the 
expression (fTUj) for A{x), we thus see that the A^s are discrete (since there are only discrete 
zeros in the Bessel functions), they must be negative and must also satisfy Eq. (fTK|) . How- 
ever, even though A(0) 7^ it follows from Eq. (jl4j) that -B(O) = 0. Therefore the sum in 
Eq. (HHD includes = with Aq = 0. 
The solution to © is thus 

00 

iV(x, t) = Co B{x) + J2 C„e^"*5(x + A„). (19) 

n=l 

Here we separated the first term explicitly in order to emphasize that it is time independent. 
Since the A^s are negative for n 7^ it follows that after a sufficiently long time the first 
term dominates, 

N{x,t) Co B{x) for t 00. (20) 

Please note that Co < 0. In order to use Eq. (fT^ we need to determine the constants C„. 
This problem will be discussed in the next section. 



III. ORTHOGONALITY 



The solution ()19|) is incomplete as it stands since we need to determine the expansion 
coefficients. This problem is related to orthogonality of the functions entering this solution. 
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It turns out that the B{x + A„)'s are not orthogonal j33|, which is connected to the fact 
that the i?— function satisfies the third order differential equation (jS)). However, we shall 
now use the fact that A{x) satisfies the second order differential equation (jl3p . so we can 
use standard orthogonality consideration as far as A{x) is concerned. Therefore instead of 
N{x,t) we consider the function 

oo 

M(x,t) = 5]C„e^"*A(x + A„), (21) 

n=0 

where the constants Cn are the same as in Eq. (jHl). The connection between N{x,t) and 
M{x, t) is given by 

dl M{x, t) = -N{x, t) (22) 

as a consequence of Eq. ((7j) . 

We can now show that the functions A{x + Xn) form an orthogonal set of functions for 
A„ 7^ 0. From Eq. (fTSj) we have 

A"{x + Xn) = ix + Xn)A{x + Xn). (23) 

This second order equation allows us to use standard methods, in contrast to the original 
third order equation for N{x,t). Thus we obtain 

roc roo 

/ dx [A{x + Xm)A"{x + Xn) - A{x + Xn)A"{x + Xm)] = (A„-Am) / dxA{x+Xn)A{x+Xm). 

JO JO 

(24) 

Following standard procedures we now shift the two differentiations in the first term in the 
integrand on the left hand side by two partial differentiations. In this way the first term 
cancels the second, except for contributions from the boundaries. For x = oo there are 
no contributions, since the Bessel function vanishes exponentially. Taking into account the 
contributions from the lower limit x = we have 

/■oo 

(Xn-Xm) dxA{x + X„)A{x + Xm) = -A{X^)A'{Xn) + A'{Xm)A{Xn). (25) 
Jo 

However, since A„ are the zeros of the function A ioi n ^ 0, we have A{Xm) = ^(A„) = 0. 
It is important to notice that neither A{0) nor A'{0) vanish. For n ^ the quantity on the 
right hand side of Eq. (j25j) vanishes, and hence 

roc 

{Xn — Xm) / dxA{x + Xn)A{x + X-m) = for ?T, and m ^ 0. (26) 

JO 



Thus the set {A{x + Xn)} consists of orthogonal functions, except for Aq = 0. It then follows 
that the expansion coefficients in Eq. ()21|) are determined through the equation 3^ 



c --L 



m 



dx M{x, 0)A{x + \m) - Co dxA{x + \m)A{x) 



, m^O, (27) 



where 



roc 

Im= dx A{x + \raf. (28) 

J 

In the next section we shall show the Im can be expressed in terms of yl'(A„), which in turn 
can be expressed in terms of Bessel functions [28] . 
Taking = we obtain from Eq. ()25|) 

r dxA{x + \n)A{x) = (29) 

Jo An 

which we shall use in the next section. 

In the expression ()27|) for C„ the function M(a;, 0) enters. However, the relevant initial 
function is N{x, 0). We can reexpress Cm in terms of N{x, 0) by solving Eq. ()22p for M in 
terms of A^. Noticing that the Greens function in one dimension is 6{x)x^ we easily obtain 
from (j22D 

M(x, 0) = a + /?x - r dx' {x - x')N{x', 0), (30) 
Jo 

where a and (3 are integration constants. It follows from the definition, Eq. (PT|) . of M that 
M(0, 0) = CoA(O), since A(A„) = Q forn ^ 0. Thus a = CoA{0). 
To find the constant (3 consider 

oo 

d,M{x,0)U=o = P=T. CnA'iXn). (31) 

n=0 

Now from the solution (jH)) we have by use of ((Zj) 

POO roo °^ ^oo 

/ dxN{x,0) = J2Cj dxB{x + \n) = -Y.Cn dxdlA{x + X^) = Y.^nA' {Xn) ■ 

•'^ n=0 •'^ n=0 •'^ 

(32) 

Therefore 

roo 

(3= dx N{x,0). (33) 
We now insert these results in Eq. fl30|) which becomes 

rx roo 

M{x,0) = - dx {x-x') N{x',0)+x dx' N{x' ,0) + CoA{0). (34) 
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To proceed we need also to determine Cq. We have by use of Eqs. (jH)) and ((Tj) 

POO POO POO 

/ xA^(x, 0) = V C„ / dx xB{x + Xn) = -^Cn dxxdlA{x + \n) = -Co A{Q), 

•^0 n=0 n=0 ''^ 

(35) 



so 

1 /-"^ TT 



Thus Co is determined in terms of the initial data for A^. 
We can now insert the result (jHU)) in Eq. (j2D), 

M{x,0)= dx' {x-x')N{x',0). (37) 

Jx 

This finally gives an expression for the expansion coefficients in terms of the initial value for 

N, 

]^ roc r POO "I 

C„ = — / dx A{x + Xm) / dx' {x - x')N{x', 0) - Co A{x) , (38) 

where we inserted the solution ()37|) of M in terms of in the expression (jTTj) for Cm- In 
Fig. 121 we show an example of the evolution of N{x, t) with the interesting feature of the 
presence of a secondary peak. Interchanging the x and x' integrations this expression can 
be rewritten in a form which is more suitable for insertions of data for N{x, 0), namely 
1 



Cm — -j— 



dx N{x, 0) / dx' {x' — x)A{x' + A^) — Cq dx A{x + Xm)A{x) 
Jo Jo 



oo 



]^ roo r rx 

— dx N{x, 0) <^ A'{x + Xm) - A'{Xm) - (x + Xm) / dx' A{x' + Xm) 

Im. JO I JO 

dx' A{x' + Xm)A{x') \ . (39) 



X 



A{0) Jo 

To obtain the last form we used Eq. (jZj) and inserted Cq. Thus, to sum up the solution of 
Eq. (jH^ is given by Eq. (|7i|l . with B given by Eqs. (fT^ and (fTIH) . and with C„ given by 
Eq. ()39|) . It should be emphasized that in Eq. (jH^ the initial function A^(x, 0) enters only 
in the first integral, whereas the other integrals involving the function A can be determined 
numerically to any desired accuracy. 

IV. ON THE EXPANSION COEFFICIENTS 

In this section we shall examine the expansion coefficients Cm given by Eqs. ()38|) and 
(jnni)- We start by considering the normalization integrals Im defined by The differential 
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2 4 6 8 

X 



FIG. 2: An example showing the time development of N{x, t) {x and t are dimensionless) with a 
secondary peak, which ultimately disappears. The specific values of the expansion parameters Cn 
are (Co, Ci, C2 . . .) = (-7, -.2, -.2, .05, -.035) . 

equation ^ for A{x), A" = xA, can be rewritten as 2A'A" - 2xAA' - A^ = -A^, i. e. 

(A'ixf - xAixf) = -Aixf, (40) 
ax ^ ^ 

so that the normalization integral becomes 

= / dxA{xf = A'(A^)^ (41) 



At7 

where we used that A{Xm) = 0. 

We can now express in terms of Bessel functions by differentiating A( x) g iven in Eq. 
()10|) and using standard functional relations for the relevant Bessel functions • Using 

A'{x) already evaluated in Eq. (|TT|l we have 

Im = ^Xl [J_|(2|A^|3/V3) - J|(2|A^|3/V3)]^. (42) 

If we use the well known asymptotic expansions of the Bessel functions as well as the 
asymptotic eigenvalue (fTTj) it is easy to see that Eq. becomes 



3\l/3 / 3x1/3 
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We have compared ()43p with the exact expression ()42p even for small m. We find that for 
m = the analytic expression (jl^ gives 4.85285, whereas the approximate result PHj) gives 
4.785. For m = 4 the corresponding numbers are 8.85737 and 8.8538. For m > 5 one gets 
the first two decimals right. 

The first expression can now be simplified by use of Eq. , 

Cm= ,0 / dxN{x,0) dx'{x'-x)A{x' + Xm)- , J dxxN{x,0), m^O. 

(44) 

From the asymptotic statements ()17|) and ()43|) we see that the last term on the right hand 
side of fl44|) behaves like m~^/^. However, in the expansion for M this is to be multiplied by 
the oscillating function A{x + A„). In any case it should be noticed that if t 7^ there is no 
convergence problem because the exponential damping factors e^™* give rapid convergence. 
Convergence problems in the form of very slow convergence may arise at the initial time, 
where the expansion coefficients are determined in terms of the initial data. 

V. MEAN VALUES 

We shall now compute the mean values x and x'^, which turn out to be given by rather 
simple expressions for large times. We have 

/•oo / rOG roo I r-oo 

x= / dxxN{x,t) / / dx N{x,t) , and = / dx x'^N{x,t) / / dx N{x,t) . (45) 
Jo I Jo Jo I Jo 

By the methods already used in Eqs. (jH^ and (jH^ we obtain 

X = -A{0) j (^A'(O) + (1/Co) £ C„e^"*A'(A„) j . (46) 
We also have by use of Eq. ((Zj) and two partial integrations 

/•oo ^ roo 

/ dx x^N{x, t) = -2 V C„e^"* / dx A{x + A„). (47) 
Jo Jo 

In the limit t — > 00 we obtain the simple results 



32/3 r(4/3) , 

X — — ~ 1 

as well as 



. . X - ^.3717, (48) 
r(2/3) ' ^ ^ 



x2 ^ 2.5758. (49) 
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Here we used that dx A{x) ~ 1.0472. The relative dispersion therefore becomes 



D'^ x^ — x"^ , ^ 

— = — ^ 0.36897. (50) 



X X 

for t —>■ oo. The dispersion is therefore smaller than x, D ^ 0.60743 x 



VI. TRANSFORMATION OF THE MODEL WITH FINITE COAGULATION TO 
A RICCATI EQUATION 

In this section we shall discuss the model, Eq. 0, with constant diffusion and finite size- 
independent fragmentation and coa gula tion kernels, / > and f3 > 0. The basic equation 



for the steady state solution is thus j30 1 







D diN{x) - fxN{x) + 2/ / dx'N{x') + (3/2 j dx'N{x')N{x - x') - 2(3N{x) = 0, (51) 
where 



(3 = (3/2 dxN{x). (52) 
Jo 

By integrating Eq. we easily get 

^iV'(O) = Tdx xN{x) - 2b' > 0, (53) 

where b = (3/ f represents a new length scale compared to the one given by the diffusion to 
fragmentation ratio, Xq = {D/ fY^^, as previously introduced. 

Eq. flSlj) is an integro-differential equation. In general we would like to solve this equation 
with the boundary conditions that A^(0) = and N{x) for x — > oo. Furthermore, we 
need of course also to invoke the condition that N{x) is positive. This is not automatically 
guaranteed from Eq. (j^Tj) . We shall see that if diffusion is absent, it is not possible to achieve 
A'^(O) = 0. The main result of this section is that in Fourier space the integro-differential 
equation (j5ip can be mapped to a Riccati equation, or alternatively, a linear second order 
differential equation. If we can achieve A^(0) = we need A^'(O) > from positivity of N{x), 
and hence the right hand side of Eq. must be positive. Since Eq. is a consequence 
of a solution for the function N{x) should satisfy automatically. 

In the presence of diffusion, we can make Eq. dimensionless by rescaling the cluster 
size X ^ x/xo as in section II: 

dlN{x) - xN{x) + 2 / dx'N{x') + ^— / dx'N{x')N{x - x') N{x) = 0, (54) 



2/Xo ^0 Xq 
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We now take the Fourier transform 



N{x) = du e'""Ar(cj) = e~^" / du e'""Ar(cj), (55) 



where the convergence factor is needed in the following. Alternatively we can consider the 
contour to be shifted slightly to the upper half plane. We assume that the function Niu) is 
continuous as the real axis is approached. 

We can now translate Eq. (j54|) to Fourier space. For example, we use 

/ dx' N{x') = / duj N{uj) / e^"" = i / duj e'^\ (56) 

Jx J—oo Jx J—oo+ie CO 

Here the convergence factor is of course important. Putting everything together we obtain 
from Eq. 

/"^" du e'''^0{uj)N{uj) = ^ / dudu' ^^"^^^^^f'^ (g-- _ e^'^) , (57) 

J-oo+ie 2/Xo J-oo+it U! — LU' ^ ' 

where the operator O is given by 

0{uj) = -uj^ + I—. (58) 

uj xq duj 

From this we obtain by an integration over x from to oo with the weight factor exp(— irx) 

. r- «^ = ^ F(r?, (59) 

J-oo+ia UJ — T 2jXq 

where F is the Hilbert transform 

Fir) = r" m 

J—oo+it T — UJ 

Eq. ()59|) can be reformulated as an equation for F. To do this, we rewrite u"^ = (u;^ — r^) 
to obtain 3^ 

/ du ^^^I^M. = [ duj(uj + t)N(uj) + t^F(t) = -tN'(0) + tN(0) + t^F(t). (61) 

J UJ — T J 

We want to solve the problem with the boundary condition A^(0) = 0, so in the following 
the linear term on the right hand side will be left out. We also have 

duj—^^ = -f + - Fir). 62 

U![U! — T) pr T 
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Collecting these results Eq. (fS^ becomes 



P dFir) (2 . 2ih\ ^, , ^ /3 

^ -F(r)2 = + _ + ^1-2 + — ^(7-) _ ]v'(0) - 2i-5-. (63) 



This is a standard equation of the Riccati type. It should be noticed that if there is no 
diffusion the Riccati equation simplifies to 

^F{Tf = -^£ill + 2tb) Fir) - 2i4- (No diffusion), (64) 

// ar \t j jjT 

where r now has the dimension of length^^. 

It is well known that a Riccati equation can be transformed to a linear second order 

differential equation by the substitution 

^ F{t) = j- log «(r), (65) 



2fxo dr 
where u satisfies 

""(^' - (7 ^ ''3 (^^''°' 1^) ^ '''' 

The original Fourier transform (j55|) was given in terms of the function A^, which in turn 
is related to the function F by the Hilbert transform ()60|) . However, it is well known that 
the Fourier transform of a Hilbert transform has a remarkably simple property, which in our 
case (since a; > 0) leads to 

N{x) = — / dTe'^''F{T). (67) 

This equation can be derived by noticing that since the u;— integration is shifted to slightly 
above the real axis, we have 

F(r) = Fp(r) + z7riV(r), Fp(r) = P T dco (68) 

J-00 T — Ul 

where P means the principal value. Further we have 

duo Niuo) e'"" = — duo Niuj) e'"" P dy = — / dre'^'' Fpir), (69) 

zvr J J~oo y iTC J 

where we used that for x > 

/oo piyos 
dy = i7T, (70) 
y 

and where we took y = t — uo. Using Eq. ()68j) we obtain Eq. (jHTj) . 

From Eq. ()67|1 we see that the original integro-differential equation ()54|1 is solved if we 
can solve the Riccati equation (j63|) . or, alternatively, the linear second order differential 
equation 
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VII. THE CASE OF NO MERGING 



Before we discuss the Riccati equation it is quite instructive to consider the case of no 
merging, d = 0. In our previous work we have solved this problem in terms of the Airy 
function [l^ ll^ with the result 

/oo 
/^+*^^ (71) 
-oo 

This function solves the equation 

0{yj)N{yj) = (72) 

subject to the boundary condition that N{x) ^ for x — ^ oo, 
With the present method we obtain from Eq. (jUHj) with j3 = 



\- IT 



F{t) - N'{0) - 2i- = 0, (73) 



dr \T J T 
where the constant c is given by c = = /q°^ dxN{x). The solution is given by 

dr 



Fir) = e^^'/3 



K-J^e-'^'/' (^iV'(0) + 2z^ 



(74) 



Here i^' is a constant of integration. Comparing to the correct result (j71|) by use of the 
Fourier transform (j67|) in terms of F, we see that the second term in the solution (|74|1 
apparently gives a deviation from ()71|). 

The resolution of this paradox is that the last term in Eq. (ff^ does not give any 
contributions to the Fourier integral. This can be seen by noticing that the indefinite integral 
in Eq. (j74j) can be expressed in terms of the incomplete F function, and when multiplied by 
the prefactor e*"^^/^ it is analytic in the upper half plane. When we plug the second term in 
(I74|) into the Fourier integral ()67|) we can close the integration by a large circle in the upper 
half plane. This circle contributes nothing, since the factor exp(irx) gives an exponential 
damping on the circle, and since the behavior of the second term (including the prefactor) 
in fj74p is —iN'{0)/T'^ + 0(l/r^) for large r, as can be seen by a partial integration. This 
behavior does not compete with the exponential damping. Therefore 

_ e^-V3 J J ^~^rV^ (^jv'(O) + 2z^) (75) 
is a null function in the Fourier integral (|67p. 
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The moral of this story is that although the function N{x) is given by the two integrals 
()55|) and (jUTjl . this does not mean that the two functions and F are proportional. They 
can differ by a non-trivial null function. As the example given above shows, we must in 
general expect the occurrence of such non-trivial null functions. 

Another point which should be mentioned is why the first term in the solution (|74j) cannot 
be treated like the second term: In the second term the factors exp(±ir^/3) cancel out, as 
one can see by repeated partial integrations of the integral in the solution (f?^ . For example, 
for small r one has 

e"'l^ [ = -- - -r' e^^'^' f dr e^^/^^ = -- - -r' + ... . (76) 

J 2 2 J 2 2 ^ ^ 

A similar argument does not apply to the first factor in (f7^ which contains exp(ir'^/3) and 
does not give convergence on the large circle: This factor overwhelms the damping from 
exp(zra;). Actually the contour can only be be closed in the upper half plane by two straight 
lines in the directions it/ 6 and 5tt/6. The closed contour running along the real axis from 
— oo to +00, a part of a circle at infinity going from the angle to 7r/6 (gives no contribution 
due to a damping factor exp(— sin(30)), < < 7r/6)), runs back to the origin along 
a straight line, and moves to the left along a line making the angle bn/Q, and is closed 
by a circle going towards the real axis (gives no contribution). On the straight lines the 
factor exp(ir'^/3) becomes strongly damped and behaves like exp(— |r|^/3). Therefore the 
oscillating integral (f7T| can be replaced by the strongly damped integral 

/■oo „ 

N{x)a=o = const / dr e-^V^-rxsm^/e sin(rx cos 7r/6). (77) 
Jo 

VIII. THE CASE OF NO DIFFUSION 

We shall now consider the case without diffusion, where the Riccati equation simplified 
to Eq. with the corresponding second order equation 

u"{t) - + 2z6) n'(r) + ^n(r) = 0. (78) 

This equation can be solved in terms of Bessel functions, 

u{t) = kr^^ e'^^ (J_|(-6r) + CAr_|(-6r)) , (79) 
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where k and C are constants. These Bessel functions can be expressed in terms of trigono- 
metric functions, 



K'-^)' K'-^j 

(80) 

For the function F we then get 

^ 2/ (1 + &r(C - z)) cos(br) - (C - &r(l + tC)) sin(&r) 
^ ^ /3 {l + CbT)cos{br) + {bT-C)sm{bT) ' ^ ^ 

where C is a constant of integration. Since F is given by a Hilbert transform, we see that 
for r ^ oo 

F{r) ^ - / dujN{uj) = - N{0), (82) 

T J T 

with correction terms starting with —iN'{0)/T'^. Thus, if A^(0) = 0, we have that F behaves 
hke Going back to the solution (jST|) . we see that it is only possible to let F go like 

1/r, which fixes the integration constant to be C = i. Then 

We can now perform the Fourier transform ()67p by means of Cauchy's theorem, and we 
obtain 

N{x) = ^ e-^/^ (84) 

It is easily verified that this function satisfies the original integro-differential equation 
without the diffusion term. The solution is identical to the one found in j^l , by noticing that 
in the case of no diffusion, b can be expressed in terms of the total "mass" M = /q°° dxxN{x) 



from Eq. as b = y/3/(2/)M. Clearly, this solution does not satisfy A^(0) = 0. Instead 
there is a piling up at x = given by the ratio of fragmentation versus merging. The more 
fragmentation one has, the more piling up comes about for small x. Also, the more merging 
one has, the less piling up at small x. These results are of course physically reasonable. 



IX. THE EFFECT OF DIFFUSION 



If diffusion is included, it will have the profound effect that it is possible to enforce the 
boundary condition A^(0) = 0. We then have to consider the full Riccati equation 
Unfortunately we have not been able to solve this equation exactly, in contrast to the no 
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diffusion case considered in the last section. Therefore we need to find an approximate 
solution valid in different r— ranges. In this connection it should be pointed out that if F{t) 
is (i) analytic in the upper half plane, and (ii) grows less than exponential on a large semi- 
circle in the upper half plane, then there is only the trivial solution N{x) = 0. This simply 
follows by the use of Eq. (j67|) by closing the contour in the upper half plane at no price, 
since there is an exponential damping from the factor exp(ixr) due to (ii). Application of 
Cauchy's theorem then gives the trivial result, since no singularities are included inside this 
contour, due to (i). Therefore, to have a non-trivial approximate solution at least one of the 
two conditions (i) and (ii) must be violated. 

To see that the boundary condition A^(0) = is possible when diffusion is included, we 
notice that for large r the linear second order equation has the solution 

u{t) ^ const. ( 1 + N'{0) + 0(l/r^) ) , i.e., logM(r) ^ const. +-^^ N'{0)+O{1/t'^) 

(85) 

which comes as a cancellation of the diffusive terms —ir'^u'^r) and +2;^ A^'(O). Using Eq. 
(jnSl) we have 

It should be noticed that in the case of no diffusion the asymptotic behavior of u{t) is 
different, u{t) ^ const. + O(logr). 

We shall now show that if one has an approximate pole behavior of F{t) this represents 
a saddle point in the Fourier integral ()67|) . This is important, since we cannot directly use 
Cauchy's theorem to perform the integral if the solution is only approximate and dominates 
only near a pole tq. From the Fourier integral we see that there is a stationary phase (or 
saddle point) when 

ixT - log(r - To) (87) 

is stationary, which happens for 

1 1 
ix = 0, i.e. r = To H . (88) 

T — Tq IX 

Also, the second derivative of the phase — is negative, so the saddle point is stable. The 
value of N{x) in this saddle point is 

N{x) ^ const. 5^e^"« ^ (89) 
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where we should sum over all relevant tq's subject to the condition that N{x) is finite and 
positive. 

We shall now investigate the possibility that there exists a pole for small values of r. To 
this end we shall use that from Eq. (j65j) a pole in F{t) can show up as a zero in the function 
u{t), which satisfies the linear second order equation . Let us tentatively assume that 
this pole occurs for small values of r, so that we can expand 

u{t) ^ 1 + CiT + C2r^ + (90) 

where the first constant is taken to be 1, since the overall scale of u{t) is irrelevant for F{t). 
From ()66|) we then obtain 

u{t) ^l + t—T + N'{Oy + ... . (91) 
xo 2xof 

There is indeed a pole u{to) ~ 0, 

= /3/ifxl)N'{0) (V(V^o)^ + 2/?/(/xo)iV'(0) - b/xo) , (92) 

where we took the root which leads to a finite result for N{x) by use of Eq. ()89|) . Using Eq. 
fISHj) this can be written 



i 2 < X > ^ \ 

"°= <x>-6/a:o W " T ^ ^ 

It is clear that this result is only valid provided |ro| is small, since otherwise more terms 
should be included in the expansion (jHOl) 

Now the question naturally appears as to whether there are other saddle points relevant 
when T is large? In general, the saddle points are in the complex r— plane, as we saw in the 
pole case discussed above, where tq is imaginary. In this connection we remind the reader 
that in the case of no merging, we obtained as a solution an oscillating integral (jTlj) . which 
could be turned into an integral (f77j) which is damped at large r by deforming the contour. 
On the deformed contour the Fourier transform is thus small for large r. Motivated by these 
considerations we ask if there exists an approximate solution of the Riccati equation ()63p 
where F is small. In this case we can ignore the quadratic term, and the equation reduces 
to 

_ f!£M + f ! + + ^J^] f (,) _ A,.(0) - 2,1 « 0, (94) 

dr \T Xo J pr 
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We see that except for the term ^F(r) the equation is the same as the equation with 
no merging (fTSj). and the solution is the same as Eq. (fTij) . except for a contribution 
exp{2ib/xoT). Again there is a null function with respect to the Fourier integral ()fj7|l . similar 
to the one that occurs in (ff^ . Ignoring the null function we have 

F(r) ^ const. e^^'+^^. (95) 

Along the contour forming angles 7r/6 and 57r/6 with the real axis this function is small for 
large values of r. The expression for N{x) is then 

N{x) ^ const, r c/rr^e^^'^^^^'''^. (96) 

J —oo 

For large x this leads to a saddle point in the Fourier transform ()(j7|l well known from the 
theory of Airy functions. It occurs at r = i^/x and is stable, 

N(x) ~ const, e~3^ . (97) 

To conclude this discussion, we see that there are in principle two saddle points (at least, 
there may be others that we have overlooked). The one given by the slope ()93p always 
dominates, because it decays exponentially in contrast to (P7|) . However, there may be a 
transitory region for large, but not too large x, where the behavior is given by (jHTj) . but for 
very large x the behavior is always dominated by the slope ()93j) . If the slope ()93j) increases 
to around one, the expression (|93|) is no longer expected to be valid, since more terms would 
be needed in the expansion (|^. 



X. NUMERICAL RESULTS 



The most direct way of validating the theoretical predictions presented in the preceding 
section is to integrate Eq. (j3T|) numerically. Inevitably, a numerical approach implies that 
the domain of the function N{x) will be restricted to a finite set of natural numbers X = 
{0, 1, ■ ■ ■ , L}. Here, we must ensure that L is chosen sufficiently large to suppress finite size 
effects and that X is sufficiently "dense" to give a proper representation of the continuous 
function A^. The second requirement is satisfied by rescaling the diffusion to fragmentation 
ratio so that Xq^ = (^j^ ^ ^1. Consequently, L must be chosen such that L/xq ^ 1. 
As we have seen, the other relevant length scale is set by the coagulation to fragmentation 
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ratio b = jj N{x). Since diffusion processes will always be dominating for small values 
of X, this length scale only plays a part in setting the upper limit, i. e. L/b ^ 1. The most 
important concern in choosing L comes from the problem of ensuring the mass conservation 
of Eq. (j51|) dt dxN{x,t)x = 0. The proper discretized version of the fragmentation and 
coagulation operators that entails mass conservation in the finite domain X reads: 

[^*^]frag [dtN]f,^^ = -f{x-l)N{x) + 2f Yl (98) 

x'=x+l 

and 

[dtNjcoag ^ [mhoag = (3/2 J2 N{x')N{x - x') - f3N{x) J2 N{S:'), (99) 

x'=0 x'=0 

where N represents the discretized function. The discrete Laplacian 



dlN{x) A^N{x) 



N{1) - 2iV(0) for X = 

N{x + 1) + N{x - 1) - 2N{x) for < i < L, (100) 

N{L - 1) - 2N{L) iov x = L 



however will always lead to a net loss of the total mass, since: 

J2xA'^N{x) = -{L + 1)N{L). (101) 

x=0 

Here, we take the simplest approach of setting L so that the relative loss of mass will be small 
in each time step, i.e. LN{L) <^ Y.^=o^N(x). Since for all parameters of b/xo we expect 
the distributions to display an exponential decay (at least), the mass loss will automatically 
be small when L/xq ^ 1 and L/b ^ 1. In the following we present results using xq ~ 50 
and L = 1000 for values of b/xo <2 and L = 10000 for b/xo > 2. The discretized equation; 

dtN{x) = A'N{x) + [dtN]i,^g + [dtN]coag (102) 

is solved using fourth order Runge-Kutta integration scheme with variable time-step Isij ]. 
Convergence is obtained by requiring that the relative change AAf of the total counts, JV = 
J2x=o^i^) during a time interval of At = 1 satisfy AJ\f /J\f < 10"^. For all calculations the 
time-integrated mass loss is less than 0.1% compared to the mass of the initial distribution. 

In Fig. Elwe show N{x) = N{x) as function of x = x/xq for different values of b/x^. For 
small values of b/xQ, one indeed observes a cross-over from the functional form Eq. (j97j) at 
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FIG. 3: The numerical steady-state distributions N{x) as function of x (both quantities dimen- 
sionless) for different coagulation-to- fragmentation ratios b/xQ. 

intermediate values of x to a pure exponential for large values of x. At larger values of 6/a;o, 
only a single exponential form is observed for x > 1. 

In Fig. |3]we show the slope IitlItq) of the exponential tail of the distributions as function 
of b/xQ. The full curve is expression flUH|) . whereas the square diamonds are the slopes 
obtained from an exponential fit to the tail of the numerical solutions for different values of 
b/xQ. For Q'(ro) < 1, which corresponds to values where the coagulation to fragmentation 
ratio is the dominant length scale, b > Xq, we find perfect agreement between the theoretical 
predictions and the numerical results. 

We shall now use a simple ID lattice model as an alternative approach to put the afore- 
mentioned theory to a test. We pick a finite lattice of length L and impose periodic boundary 
conditions. On the lattice we initially distribute a set of Af clusters of sizes {ii}iLi and with 
a total mass chosen to fit the lattice length, J2i = L. We define the cluster boundaries by 
a new set {xi}^^ such that ii = Xi — Xi^i (with xq=xn — L). The diffusion process of Eq. 
()51|) corresponds in the lattice model to a random walk of the boundaries and therefore in 
each time step we update the boundary positions according to a;i(t + 1) = Xiit) ± 1. Here 
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FIG. 4: The slope /m(ro) (dimensionless) of the exponential tail of N(x) as function of b/xQ. The 
solid line is the theoretical prediction, Eq. (jHSI) and the square diamonds are the numerical results. 
One observes a perfect agreement for all /m(ro) < 1. 

the lattice spacing and time steps, as well as the diffusion constant, are chosen to be unity. 
A fragmentation event is simulated by introducing a new random walker or boundary x on 
the lattice. The position of x is picked uniformly among the L — M available lattice sites. 
Effectively, we therefore get a constant fragmentation rate given by / = n/ where n is the 
number of new boundaries introduced in each time step. Finally, the coagulation process is 
similar to the removal of boundaries. In each time step a random walker is removed with a 
probability 7. Note that 7 is equivalent to the rate (3 = (3M introduced in (j52|) . The equiva- 
lence follows by comparing the number of events where a cluster of size x coagulates with a 
cluster of size y. In the lattice model, the number of events in a unit time are J^l^^^^^^^, 
which should be compared with the corresponding number in Eq. ()51|) . and thus we get 
7/A/' = (3. In order to minimize the effect of the underlying lattice it is important to keep 
the total mass or lattice length L of the system much larger than any other scale in the 
system, i.e. L ^ {D/fY^^ and L ^ 13/f. Fig. Elpresents estimates of (j93|) based on the 
lattice model with a lattice size L = 10''. We collect the simulation data in a histogram 
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FIG. 5: Lattice simulation data of Imro as function of 6 = (3/ f ■ The line represents the analytical 
expression of Imro- 

with a unit bin-size and then fit the exponent of the taiL The maximum count in any bin 
is for reasonable amounts of computer time hmited above by 10^, thus we have a bound on 
the maximum exponent that it is possible to estimate and we have an explanation why the 
lattice simulations provide poor estimates in Fig. El for low values of h. 

Compared to the direct numerical integration there are both advantages and disadvan- 
tages. First of all the lattice model is extremely simple from a computationally point of 
view and does by construction conserve the total mass. Moreover the model reproduces to 
high accuracy, using little computer time, the theoretical predictions for large values of 6/xo, 
which is in contrast to the larger computational power needed in the direct numerical inte- 
gration. Note that in the lattice model there is a weak size correlation between neighboring 
clusters. The diffusion makes some clusters large on the cost of the surrounding ones. The 
correlation is not present in the mean field equation ()5H1 and it may be avoided on the lattice 
by shuffling in each time step the clusters. By implementing this shuffling in the numerical 
routine, it however turns out, that the correlation has negligible or no effect on the cluster 
size distribution. 
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XI. CONCLUSION 



In this paper, we have introduced a model which includes the three fundamental physical 
processes: diffusion, fragmentation and coagulation. The model is formulated as a dynamical 
equation in terms of the distribution N{x, t) of fragment sizes x. The main results of the 
paper are: In the case of no coagulation term, we obtain an exact solution for the distribution 



non-linear equation can be mapped exactly onto a Riccati equation. From solution of this 
equation we obtain that the distribution now turns into a pure exponential for large x. 
When the coagulation process is small as compared to the fragmentation process, we identify 
directly that the distribution N{x) behaves as the Bessel function for small x after which it 
crosses over to an exponential at a specific value of x. 

We believe that our proposed model is relevant in many physical situations, such as for 
solutions of macro molecules like polymers, proteins and micelles. In fact, from measured 
distributions of fragment sizes we suggest that one might be able to identify how important 
the coagulation process is compared to the fragmentation process. We are in the process of 
collecting experimental data for such investigations. 
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